import numpy as np
import pandas as pd
import scipy.integrate
import scipy.interpolate
import scipy.optimize
import biocircuits
import bokeh.io
bokeh.io.output_notebook()
import panel as pn
pn.extension()
dark_green = "#31a354"
light_green = "#74c476"
dark_orange = "#e6550d"
light_orange = "#fd8d3c"
def style(p, autohide=False):
p.title.text_font="Helvetica"
p.title.text_font_size="16px"
p.title.align="center"
p.xaxis.axis_label_text_font="Helvetica"
p.yaxis.axis_label_text_font="Helvetica"
p.xaxis.axis_label_text_font_size="13px"
p.yaxis.axis_label_text_font_size="13px"
p.background_fill_alpha = 0
if autohide: p.toolbar.autohide=True
return p
There is a small typo in what's below, in part (b), df/da should be scaled by a $\kappa$ term, but this does not effect the analysis in part (b), and corresponding changes have been made to the code in the later parts.

We can show that the fixed point is stable by plotting a phase portrait and if the derivative has a negative slope during the zero-crossing, it means that when [$A$] is below [$A_\mathrm{st}$], the derivative will be positive and [$A$] will increase until it reaches $A_\mathrm{st}$, and when [$A$] is above [$A_\mathrm{st}$], the derivative will be negative and push the concentration to $A_\mathrm{st}$
def analytical_steady_state(alpha, beta, gamma, kappa):
_term1 = gamma*(alpha-1)
_term2 = kappa*(beta+gamma)
_term3 = np.sqrt( 4*alpha*gamma*(kappa+gamma) +
(gamma*(alpha-1)+kappa*(beta+alpha))**2 )
numerator = _term1 + _term2 + _term3
denominator = 2*gamma*(gamma+kappa)
A_ss = numerator / denominator
R_ss = gamma * A_ss
return A_ss, R_ss
# function for integrations
def derivs(AR, t, alpha, beta, gamma, kappa, n):
A, R = AR
f = (kappa*A)**n / (1 + (kappa*A)**n + R**n)
deriv_A = alpha + beta * f - gamma*A
deriv_R = alpha + beta * f - R
return np.array([deriv_A, deriv_R])
# function for plotting
def find_derivs(A, R, alpha, beta, gamma, kappa):
f = kappa*A / (1 + kappa*A + R)
deriv_A = alpha + beta * f - gamma*A
deriv_R = alpha + beta * f - R
return deriv_A, deriv_R
n = 1
alpha = 1
beta = 1
gamma = 2
kappa = 1
A_ss, R_ss = analytical_steady_state(alpha, beta, gamma, kappa)
_A_range = np.linspace(0, max(A_ss, R_ss)*2, 100)
_R_range = np.linspace(0, max(A_ss, R_ss)*2, 100)
_A_deriv, _R_deriv = find_derivs(_A_range, _R_range, alpha, beta, gamma, kappa)
_A_ss_deriv, _R_ss_deriv = find_derivs(A_ss, R_ss, alpha, beta, gamma, kappa)
p = bokeh.plotting.figure(height=350, width=400, title="Phase Portrait",
x_axis_label="[A], [R]", y_axis_label="d/dt")
q = bokeh.plotting.figure(height=350, width=400)
# zero line
p.line(( min(_A_range.min(), _R_range.min()),
max(_A_range.max(), _R_range.max())), 0,
line_dash="dotdash", color="black", line_width=3, alpha=0.5 )
# .... PHASE PORTRAIT ....
f = kappa*A_ss / (1 + kappa*A_ss + R_ss)
p.line(_A_range, _A_deriv, line_color=dark_green, line_width=3)
p.line(_R_range, _R_deriv, line_color=dark_orange, line_width=3)
# steady state vertical lines
p.line(A_ss, ( min(_A_deriv.min(), _R_deriv.min()),
max(_A_deriv.max(), _R_deriv.max()) ),
line_dash="dashdot", color=light_green, line_width=3 )
p.line(R_ss, ( min(_A_deriv.min(), _R_deriv.min()),
max(_A_deriv.max(), _R_deriv.max()) ),
line_dash="dotdash", color=light_orange, line_width=3 )
# .... INTEGRATING TRAJECTORIES ....
N_SEED = 50
scale = 2*max(A_ss, R_ss)
t = np.linspace(0, 6, 500)
args = (alpha, beta, gamma, kappa, n)
ARos = np.random.random((N_SEED,2))*scale # random initial conditions
q = bokeh.plotting.figure(height=350, width=400, title="activator A", x_axis_label="time", y_axis_label="[ A ]")
r = bokeh.plotting.figure(height=350, width=400, title="repressor R", x_axis_label="time", y_axis_label="[ R ]")
for ARo in ARos:
_AR = scipy.integrate.odeint(derivs, ARo, t, args=args)
A, R = _AR.T
q.line(t, A, color=dark_green, line_width=2, line_alpha=0.6)
r.line(t, R, color=dark_orange, line_width=2, line_alpha=0.6)
# q.line(t, A_ss, color=light_green, line_width=5, line_dash="dotdash")
# r.line(t, R_ss, color=light_orange, line_width=5, line_dash="dotdash")
bokeh.io.show(bokeh.layouts.layout([[style(p), style(q), style(r)]]))
On the left, we see the negative slopes which indicate stability. For fun, we can also seed 50 initial random conditions (random in both [Ao] and [Ro]), integrate their trajectories, and watch them converge to the stable points. Nice little horizontal funnels of joy.
Do we need ultransensitivity? Here, we've shown that with n = 1 and leakage, there is only one stable fixed point, and sustained oscillations cannot occur. Ultrasensitivity would help. However, I am having trouble understanding how ultrasensitivity would manifest for single-occupancy. In the very first week, I remember discussions about how an additional cooperative binding term could result in higher order terms, but for single occupancy, in my mind all the binding sites are operating independently of each other, so how would this work?
By Des Carte's rule of sign, 3 sign changes shows either a maximum of 1 or 3 fixed points on the positive domain of a. So at most we have 3, and if we can show that we have 3, then that is the maximum. Let's plot this fifth order polynomial in a desmos plot (note that all the letters are alphabetical coefficients, and $a$ is $x$). How I went about the fixed points was first choosing the last two coefficients ($\gamma$ and $\alpha$) first, both to be less than 1 (since the first two coefficients are exponentiating $\gamma$, and scaling by $\alpha$ in the second), and then finding the first two coefficients that satisfied 1, 2, 3, then back-solving for $\kappa$ from the first coefficient, and plugging that into the second coefficient for $\beta$.
Note that the plot of the polynomial is not the derivative and will not provide information about stability, but instead takes the steady state condition as a polynomial and shows the zero-s.
The parameters are shown in the table below. For each, I will plot the polynomial to visually check the zero points. Then, we will drop a bunch of initial conditions and see how they behave.
| # of fixed points | $\alpha$ | $\beta$ | $\gamma$ | $\kappa$ |
|---|---|---|---|---|
| 1 | 0.69 | 0.59 | 0.79 | 1.21 |
| 2 | 0.12 | 0.51 | 0.49 | 1.q65 |
| 3 | 0.12 | 1.26 | 0.79 | 1.21 |
The trajectories take ~5 seconds to run.
def plot_polynomial(alpha, beta, gamma, kappa, n, title, low=0, high=10):
coeff1 = gamma*(gamma**n + kappa**n)
coeff2 = (kappa**n)*(alpha + beta) + alpha*gamma**n
coeff3 = gamma
coeff4 = alpha
x = np.linspace(low, high, 500)
y = coeff1*x**(n+1) - coeff2*x**n + coeff3*x - coeff4
p = bokeh.plotting.figure(height=250, width=400, title=title,
x_axis_label="A", y_axis_label=f"{n+1}th order poly")
p.line(x, y, line_color="#40312b", line_width=3)
p.line(x, 0, line_color="#40312b", line_width=3, line_dash="dotdash")
return style(p, autohide=True)
def fixed_point(alpha, beta, gamma, kappa, n, low=0, high=100):
def _root_function(a):
coeff1 = gamma*(gamma**n + kappa**n)
coeff2 = (kappa**n)*(alpha + beta) + alpha*gamma**n
coeff3 = gamma
coeff4 = alpha
return coeff1*a**(n+1) - coeff2*a**n + coeff3*a - coeff4
return scipy.optimize.brentq(_root_function, low, high)
n = 4
alpha = 0.69
gamma = 0.79
kappa = (2/gamma - gamma**n)**(1/n)
beta = (3-alpha*gamma**n)/(kappa**n)-alpha
fp1_i = fixed_point(alpha, beta, gamma, kappa, n)
p = plot_polynomial(alpha, beta, gamma, kappa, n, f"n = 4: 1 fixed point", low=-0.1, high=1.7)
alpha = 0.12
gamma = 0.49
kappa = (3.7/gamma - gamma**n)**(1/n)
beta = (4.72-alpha*gamma**n)/(kappa**n)-alpha
fp1_ii = fixed_point(alpha, beta, gamma, kappa, n)
fp2_ii = fixed_point(alpha, beta, gamma, kappa, n, low=1)
q = plot_polynomial(alpha, beta, gamma, kappa, n, f"n = 4: 2 fixed points", low=-0.1, high=1.3)
alpha = 0.12
gamma = 0.79
kappa = (2/gamma - gamma**n)**(1/n)
beta = (3-alpha*gamma**n)/(kappa**n)-alpha
fp1_iii = fixed_point(alpha, beta, gamma, kappa, n)
fp2_iii = fixed_point(alpha, beta, gamma, kappa, n, low=0.16, high=1)
fp3_iii = fixed_point(alpha, beta, gamma, kappa, n, low=1, high=2)
r = plot_polynomial(alpha, beta, gamma, kappa, n, f"n = 4: 3 fixed points", low=-0.1, high=1.5)
# trajectories...
plots = []
for n_fps in [1, 2, 3]:
N_SEED = 250
tmax = 20
if n_fps == 1:
fp1, fp2, fp3 = fp1_i, -np.inf, -np.inf
alpha, gamma = 0.69, 0.79
kappa = (2/gamma - gamma**n)**(1/n)
beta = (3-alpha*gamma**n)/(kappa**n)-alpha
if n_fps == 2:
fp1, fp2, fp3 = fp1_ii, fp2_ii, -np.inf
alpha, gamma = 0.12, 0.47
kappa = (3.7/gamma - gamma**n)**(1/n)
beta = (4.72-alpha*gamma**n)/(kappa**n)-alpha
tmax = 50
if n_fps == 3:
fp1, fp2, fp3 = fp1_iii, fp2_iii, fp3_iii
alpha, gamma = 0.12, 0.79
kappa = (2/gamma - gamma**n)**(1/n)
beta = (3-alpha*gamma**n)/(kappa**n)-alpha
scale = 2*max(fp1, fp2, fp3, gamma*fp1, gamma*fp2, gamma*fp3)
t = np.linspace(0, tmax, 500)
args = (alpha, beta, gamma, kappa, n)
ARos = np.random.random((N_SEED, 2))*scale # random initial conditions
_p = bokeh.plotting.figure(height=350, width=400, title="activator A",
x_axis_label="time", y_axis_label="[ A ]")
_q = bokeh.plotting.figure(height=350, width=400, title="repressor R",
x_axis_label="time", y_axis_label="[ R ]")
for ARo in ARos:
_AR = scipy.integrate.odeint(derivs, ARo, t, args=args)
A, R = _AR.T
_p.line(t, A, color=dark_green, line_width=2, line_alpha=0.3)
_q.line(t, R, color=dark_orange, line_width=2, line_alpha=0.3)
plots.append(bokeh.layouts.layout([style(_p), style(_q)]))
bokeh.io.show(bokeh.layouts.layout([[p, q, r]]))
print(f"fixed points @ A = \t {np.round(fp1_i, 2)} \t\t\t\t\t {np.round(fp1_ii, 2)}, {np.round(fp2_ii, 2)}\t\t\t\t\t {np.round(fp1_iii, 2)}, {np.round(fp2_iii, 2)}, {np.round(fp3_iii, 2)}")
pn.Row(
plots[0],
plots[1],
plots[2]
)
fixed points @ A = 1.45 0.32, 1.22 0.15, 0.75, 1.36